Intro

Main analysis: This scripts reads spatially referenced linear trends over time in biomass, metabolic index and temperature (“03_calculate_velocities.Rmd”), and fits spatial sdmTMB models to those

# Load libraries, install if needed
library(tidyverse)
#> Warning: package 'tidyr' was built under R version 4.0.5
library(tidylog)
library(RCurl)
library(devtools)
library(sdmTMB)
library(rMR)
library(patchwork)
library(raster)
library(VoCC)
library(RColorBrewer)
#> Warning: package 'RColorBrewer' was built under R version 4.0.5
library(viridis)
library(ggh4x)
library(ggridges)

# Source code for plots
source_url("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/R/functions/map-plot.R")
#> Warning: package 'sf' was built under R version 4.0.5

theme_set(theme_plot())

pal <- brewer.pal(name = "Dark2", n = 8)

Read and plot data

# Read from GH?
d <- read_csv("data/clean/velocities.csv") %>% 
  mutate(quarter_f = as.factor(quarter))
#> Rows: 87885 Columns: 20
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr  (4): id, species, life_stage, group
#> dbl (16): X, Y, quarter, mean_log_biomass, mean_temp, mean_oxy, mean_phi, me...
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> mutate: new variable 'quarter_f' (factor) with 2 unique values and 0% NA

glimpse(d)
#> Rows: 87,885
#> Columns: 21
#> $ id                  <chr> "320.016716551297_6028.58668110674_cod_adult_1", "…
#> $ X                   <dbl> 320.0167, 320.0167, 320.0167, 320.0167, 320.0167, …
#> $ Y                   <dbl> 6028.587, 6028.587, 6028.587, 6028.587, 6028.587, …
#> $ species             <chr> "cod", "cod", "cod", "cod", "dab", "dab", "dab", "…
#> $ life_stage          <chr> "adult", "adult", "juvenile", "juvenile", "adult",…
#> $ quarter             <dbl> 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1,…
#> $ mean_log_biomass    <dbl> 4.3916628, 4.5915365, 2.8446881, 2.9603249, 3.6929…
#> $ mean_temp           <dbl> 3.369871, 10.260329, 3.369871, 10.260329, 3.369871…
#> $ mean_oxy            <dbl> 8.114926, 6.209476, 8.114926, 6.209476, 8.114926, …
#> $ mean_phi            <dbl> 1.5023498, 0.8518787, 4.6064940, 2.6120242, 4.4595…
#> $ mean_log_biomass_ct <dbl> 1.545992832, 1.745866571, -0.000981832, 0.11465494…
#> $ mean_temp_ct        <dbl> -2.888372, 4.002085, -2.888372, 4.002085, -2.88837…
#> $ mean_oxy_ct         <dbl> 2.3313114, 0.4258614, 2.3313114, 0.4258614, 2.3313…
#> $ mean_phi_ct         <dbl> -1.6917880, -2.3422591, 1.4123562, -0.5821135, 1.2…
#> $ bio_vel             <dbl> -0.03357378, -0.03357378, -0.03601321, -0.03601321…
#> $ phi_vel             <dbl> 0.3286411, 0.5693023, 0.3286411, 0.5693023, 0.3286…
#> $ temp_vel            <dbl> 1.6318561, 0.7592320, 1.6318561, 0.7592320, 1.6318…
#> $ group               <chr> "cod_adult", "cod_adult", "cod_juvenile", "cod_juv…
#> $ phi_vel_sc          <dbl> 0.8134520, 0.9991260, 0.8942800, 1.1221752, 1.0096…
#> $ temp_vel_sc         <dbl> -0.23950578, -0.47555486, -0.11759968, -0.39684290…
#> $ quarter_f           <fct> 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1, 4, 1,…

# Biotic velocities
plot_map + 
  geom_raster(data = d, aes(X*1000, Y*1000, fill = bio_vel)) + 
  facet_grid(life_stage ~ species) + 
  #scale_fill_viridis() + 
  scale_fill_gradient2() + 
  geom_sf(size = 0.3)
#> Warning: Removed 2425 rows containing missing values (geom_raster).


# Phi velocities
plot_map + 
  geom_raster(data = d, aes(X*1000, Y*1000, fill = phi_vel)) + 
  facet_grid(life_stage ~ species) + 
  #scale_fill_viridis() + 
  scale_fill_gradient2() + 
  geom_sf(size = 0.3)
#> Warning: Removed 2425 rows containing missing values (geom_raster).


# Temperature velocities
plot_map + 
  geom_raster(data = d, aes(X*1000, Y*1000, fill = temp_vel)) + 
  facet_grid(life_stage ~ species) + 
  #scale_fill_viridis() + 
  scale_fill_gradient2() + 
  geom_sf(size = 0.3)
#> Warning: Removed 2425 rows containing missing values (geom_raster).


# Plot biotic vs phi and temperature velocities
ggplot(d, aes(phi_vel, bio_vel)) + 
  geom_point(alpha = 0.4) +
  facet_grid(life_stage ~ species) + 
  stat_smooth(method = "lm") +
  theme(aspect.ratio = 1) +
  NULL
#> `geom_smooth()` using formula 'y ~ x'


ggplot(d, aes(temp_vel, bio_vel)) + 
  geom_point(alpha = 0.4) +
  facet_grid(life_stage ~ species) + 
  stat_smooth(method = "lm") +
  theme(aspect.ratio = 1) +
  NULL
#> `geom_smooth()` using formula 'y ~ x'


# Only cod
head(d) %>% as.data.frame()
#>                                                 id        X        Y species
#> 1    320.016716551297_6028.58668110674_cod_adult_1 320.0167 6028.587     cod
#> 2    320.016716551297_6028.58668110674_cod_adult_4 320.0167 6028.587     cod
#> 3 320.016716551297_6028.58668110674_cod_juvenile_1 320.0167 6028.587     cod
#> 4 320.016716551297_6028.58668110674_cod_juvenile_4 320.0167 6028.587     cod
#> 5    320.016716551297_6028.58668110674_dab_adult_1 320.0167 6028.587     dab
#> 6    320.016716551297_6028.58668110674_dab_adult_4 320.0167 6028.587     dab
#>   life_stage quarter mean_log_biomass mean_temp mean_oxy  mean_phi
#> 1      adult       1         4.391663  3.369871 8.114926 1.5023498
#> 2      adult       4         4.591537 10.260329 6.209476 0.8518787
#> 3   juvenile       1         2.844688  3.369871 8.114926 4.6064940
#> 4   juvenile       4         2.960325 10.260329 6.209476 2.6120242
#> 5      adult       1         3.692946  3.369871 8.114926 4.4595684
#> 6      adult       4         3.923447 10.260329 6.209476 2.5287129
#>   mean_log_biomass_ct mean_temp_ct mean_oxy_ct mean_phi_ct     bio_vel
#> 1         1.545992832    -2.888372   2.3313114  -1.6917880 -0.03357378
#> 2         1.745866571     4.002085   0.4258614  -2.3422591 -0.03357378
#> 3        -0.000981832    -2.888372   2.3313114   1.4123562 -0.03601321
#> 4         0.114654948     4.002085   0.4258614  -0.5821135 -0.03601321
#> 5         0.847275563    -2.888372   2.3313114   1.2654306  0.21900577
#> 6         1.077776951     4.002085   0.4258614  -0.6654249  0.21900577
#>     phi_vel temp_vel        group phi_vel_sc temp_vel_sc quarter_f
#> 1 0.3286411 1.631856    cod_adult   0.813452  -0.2395058         1
#> 2 0.5693023 0.759232    cod_adult   0.999126  -0.4755549         4
#> 3 0.3286411 1.631856 cod_juvenile   0.894280  -0.1175997         1
#> 4 0.5693023 0.759232 cod_juvenile   1.122175  -0.3968429         4
#> 5 0.3286411 1.631856    dab_adult   1.009634  -0.2510603         1
#> 6 0.5693023 0.759232    dab_adult   1.335121  -0.5594185         4

d_cod <- d %>% 
  filter(group == "cod_adult")
#> filter: removed 70,544 rows (80%), 17,341 rows remaining

p1 <- plot_map + 
  geom_raster(data = d_cod, aes(X*1000, Y*1000, fill = bio_vel)) + 
  scale_fill_gradient2(name = "Biotic velocity") + 
  geom_sf(size = 0.3)

p2 <- plot_map + 
  geom_raster(data = d_cod, aes(X*1000, Y*1000, fill = temp_vel)) + 
  scale_fill_gradient2(name = "Temperature velocity") + 
  geom_sf(size = 0.3)

p3 <- plot_map + 
  geom_raster(data = d_cod, aes(X*1000, Y*1000, fill = phi_vel)) + 
  scale_fill_gradient2(name = "\u03C6 velocity") + 
  geom_sf(size = 0.3)

((p2 + p3) / (p1 + plot_spacer()))
#> Warning: Raster pixels are placed at uneven horizontal intervals and will be
#> shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven vertical intervals and will be
#> shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven horizontal intervals and will be
#> shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven vertical intervals and will be
#> shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven horizontal intervals and will be
#> shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven vertical intervals and will be
#> shifted. Consider using geom_tile() instead.

ggsave("figures/map_velos_cod.pdf", width = 17, height = 17, units = "cm", device = cairo_pdf)
#> Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted. Consider using geom_tile() instead.
#> Raster pixels are placed at uneven vertical intervals and will be shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven horizontal intervals and will be
#> shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven vertical intervals and will be
#> shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven horizontal intervals and will be
#> shifted. Consider using geom_tile() instead.
#> Warning: Raster pixels are placed at uneven vertical intervals and will be
#> shifted. Consider using geom_tile() instead.

Fit sdmTMB models of biomass velocities as a function of phi velocities

coefs <- list()
pred_velo1 <- list()
pred_velo2 <- list()

for(i in unique(d$group)){
  
  dd <- filter(d, group == i)
  
  mesh <- make_mesh(dd, c("X", "Y"), n_knots = 200)
  plot(mesh)
  
  m <- sdmTMB(
    bio_vel ~ phi_vel_sc * mean_phi_ct + mean_log_biomass_ct + quarter_f,
    mesh = mesh,
    data = dd,
    family = gaussian(), 
    spatial = "on",
    spatiotemporal = "off",
    control = sdmTMBcontrol(newton_loops = 1)
  )
  
  sanity(m) 
  
  # Save coefficients??
  coefs[[i]] <- tidy(m, conf.int = TRUE) %>% mutate(group = i)
  
  # Now predict the biotic velocity for high and low temperature with a sd change in the climate velocity
  low_mean_clim <- as.numeric(quantile(dd$mean_phi_ct, prob = 0.05))
  high_mean_clim <- as.numeric(quantile(dd$mean_phi_ct, prob = 0.95))
  #sd_clim_v <- sd(dd$phi_vel)
  mean_log_biomass_ct <- mean(dd$mean_log_biomass_ct)
  
  # Predict
  nd <- data.frame(phi_vel_sc = c(-1, 0, -1, 0),
                   mean_phi_ct = rep(c(low_mean_clim, high_mean_clim), each = 2),
                   mean_log_biomass_ct = rep(mean_log_biomass_ct, 4),
                   X = mean(dd$X),
                   Y = mean(dd$Y),
                   quarter_f = as.factor(4))
  
  nsim = 1000
  pred <- predict(m, newdata = nd, se_fit = TRUE, nsim = nsim)
  
  pred <- data.frame(
    mean_phi_ct = rep(c(low_mean_clim, high_mean_clim), each = nsim),
    est_change = c(pred[1, ], pred[3, ]),
    est_0 = c(pred[2, ], pred[4, ]))
  
  pred <- pred %>% 
    mutate(bio_trend_delta = est_change - est_0, # if bio trend is positive, bio increases with declines in phi
           mean_clim = ifelse(mean_phi_ct == low_mean_clim, 
                              "low_phi", "high_phi"),
           group = i)
    
  pred_velo1[[i]] <- pred
  
  # Now make a prediction across a range of velocities
  # Again I'm predicting for the range of the specific model
  nd3 <- data.frame(expand.grid(phi_vel_sc = seq(quantile(dd$phi_vel_sc, probs = 0.05),
                                                 quantile(dd$phi_vel_sc, probs = 0.95),
                                                 length.out = 50),
                                mean_phi_ct = c(low_mean_clim, high_mean_clim))) %>% 
    mutate(mean_log_biomass_ct = mean_log_biomass_ct,
           X = mean(dd$X),
           Y = mean(dd$Y),
           quarter_f = as.factor(4))
  
  pred3 <- predict(m, newdata = nd3, se_fit = TRUE)
  
  pred3 <- pred3 %>% 
    mutate(lwr = est - 1.96*est_se,
           upr = est + 1.96*est_se, 
           group = i)
  
  pred_velo2[[i]] <- pred3
  
  ggplot(pred3, aes(phi_vel_sc, est, color = factor(mean_phi_ct))) + geom_line()
  
}
#> filter: removed 70,544 rows (80%), 17,341 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'bio_trend_delta' (double) with 2,000 unique values and 0% NA
#>         new variable 'mean_clim' (character) with 2 unique values and 0% NA
#>         new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#>         new variable 'X' (double) with one unique value and 0% NA
#>         new variable 'Y' (double) with one unique value and 0% NA
#>         new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#>         new variable 'upr' (double) with 100 unique values and 0% NA
#>         new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 72,173 rows (82%), 15,712 rows remaining

#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'bio_trend_delta' (double) with 2,000 unique values and 0% NA
#>         new variable 'mean_clim' (character) with 2 unique values and 0% NA
#>         new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#>         new variable 'X' (double) with one unique value and 0% NA
#>         new variable 'Y' (double) with one unique value and 0% NA
#>         new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#>         new variable 'upr' (double) with 100 unique values and 0% NA
#>         new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 86,033 rows (98%), 1,852 rows remaining

#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'bio_trend_delta' (double) with 2,000 unique values and 0% NA
#>         new variable 'mean_clim' (character) with 2 unique values and 0% NA
#>         new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#>         new variable 'X' (double) with one unique value and 0% NA
#>         new variable 'Y' (double) with one unique value and 0% NA
#>         new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#>         new variable 'upr' (double) with 100 unique values and 0% NA
#>         new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 86,343 rows (98%), 1,542 rows remaining

#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'bio_trend_delta' (double) with 2,000 unique values and 0% NA
#>         new variable 'mean_clim' (character) with 2 unique values and 0% NA
#>         new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#>         new variable 'X' (double) with one unique value and 0% NA
#>         new variable 'Y' (double) with one unique value and 0% NA
#>         new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#>         new variable 'upr' (double) with 100 unique values and 0% NA
#>         new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 64,646 rows (74%), 23,239 rows remaining

#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'bio_trend_delta' (double) with 2,000 unique values and 0% NA
#>         new variable 'mean_clim' (character) with 2 unique values and 0% NA
#>         new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#>         new variable 'X' (double) with one unique value and 0% NA
#>         new variable 'Y' (double) with one unique value and 0% NA
#>         new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#>         new variable 'upr' (double) with 100 unique values and 0% NA
#>         new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 79,901 rows (91%), 7,984 rows remaining

#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'bio_trend_delta' (double) with 2,000 unique values and 0% NA
#>         new variable 'mean_clim' (character) with 2 unique values and 0% NA
#>         new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#>         new variable 'X' (double) with one unique value and 0% NA
#>         new variable 'Y' (double) with one unique value and 0% NA
#>         new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#>         new variable 'upr' (double) with 100 unique values and 0% NA
#>         new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 79,755 rows (91%), 8,130 rows remaining

#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'bio_trend_delta' (double) with 2,000 unique values and 0% NA
#>         new variable 'mean_clim' (character) with 2 unique values and 0% NA
#>         new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#>         new variable 'X' (double) with one unique value and 0% NA
#>         new variable 'Y' (double) with one unique value and 0% NA
#>         new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#>         new variable 'upr' (double) with 100 unique values and 0% NA
#>         new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 75,800 rows (86%), 12,085 rows remaining

#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'bio_trend_delta' (double) with 2,000 unique values and 0% NA
#>         new variable 'mean_clim' (character) with 2 unique values and 0% NA
#>         new variable 'group' (character) with one unique value and 0% NA
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#>         new variable 'X' (double) with one unique value and 0% NA
#>         new variable 'Y' (double) with one unique value and 0% NA
#>         new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'lwr' (double) with 100 unique values and 0% NA
#>         new variable 'upr' (double) with 100 unique values and 0% NA
#>         new variable 'group' (character) with one unique value and 0% NA


phi_coefs <- dplyr::bind_rows(coefs) %>% 
  mutate(sig = "N",
         sig = ifelse(estimate < 0 & conf.high < 0, "Y", sig),
         sig = ifelse(estimate > 0 & conf.low > 0, "Y", sig))
#> mutate: new variable 'sig' (character) with 2 unique values and 0% NA

phi_pred_delta <- dplyr::bind_rows(pred_velo1) %>% 
  separate(group, c("species", "life_stage"),
           sep = "_", remove = FALSE) %>% 
  mutate(species = str_to_title(species),
         life_stage = str_to_title(life_stage))
#> mutate: changed 16,000 values (100%) of 'species' (0 new NA)
#>         changed 16,000 values (100%) of 'life_stage' (0 new NA)

phi_pred <- dplyr::bind_rows(pred_velo2) %>% 
  separate(group, c("species", "life_stage"),
           sep = "_", remove = FALSE) %>% 
  mutate(species = str_to_title(species),
         life_stage = str_to_title(life_stage))
#> mutate: changed 800 values (100%) of 'species' (0 new NA)
#>         changed 800 values (100%) of 'life_stage' (0 new NA)
plot_map + 
  geom_raster(data = d, aes(X*1000, Y*1000, fill = bio_vel)) + 
  facet_grid(life_stage ~ species) + 
  #scale_fill_viridis() + 
  scale_fill_gradient2() + 
  geom_sf(size = 0.3)
#> Warning: Removed 2425 rows containing missing values (geom_raster).

write_csv(phi_coefs, "output/phi_velo_coefs.csv")
write_csv(phi_pred_delta, "output/pred_phi_velo_delta.csv")
write_csv(phi_pred, "output/pred_phi_velo.csv")

Plot phi velocity models

phi_pred_delta2 <- phi_pred_delta %>%
  mutate(mean_clim = ifelse(mean_clim == "high_phi", "High \u03C6", "Low \u03C6")) %>%
  mutate(id2 = paste(species, life_stage))
#> mutate: changed 16,000 values (100%) of 'mean_clim' (0 new NA)
#> mutate: new variable 'id2' (character) with 8 unique values and 0% NA
  
order <- phi_pred_delta2 %>%
  group_by(id2, mean_clim) %>% 
  summarise(mean_delta = mean(bio_trend_delta)) %>% 
  pivot_wider(names_from = mean_clim, values_from = mean_delta) %>% 
  mutate(diff = `Low φ` - `High φ`) %>% 
  arrange(desc(diff)) %>% 
  mutate(Expected = ifelse(`Low φ` < `High φ`, "Yes", "No"))
#> group_by: 2 grouping variables (id2, mean_clim)
#> summarise: now 16 rows and 3 columns, one group variable remaining (id2)
#> pivot_wider: reorganized (mean_clim, mean_delta) into (High φ, Low φ) [was 16x3, now 8x3]
#> mutate (grouped): new variable 'diff' (double) with 8 unique values and 0% NA
#> mutate (grouped): new variable 'Expected' (character) with 2 unique values and 0% NA

order
#> # A tibble: 8 × 5
#> # Groups:   id2 [8]
#>   id2               `High φ` `Low φ`    diff Expected
#>   <chr>                <dbl>   <dbl>   <dbl> <chr>   
#> 1 Dab Juvenile       0.0154   0.0929  0.0775 No      
#> 2 Dab Adult          0.0252   0.0382  0.0130 No      
#> 3 Cod Juvenile       0.0226   0.0335  0.0108 No      
#> 4 Cod Adult          0.00753 -0.0102 -0.0177 Yes     
#> 5 Flounder Adult     0.0293  -0.0118 -0.0412 Yes     
#> 6 Flounder Juvenile  0.101   -0.123  -0.224  Yes     
#> 7 Plaice Juvenile    0.141   -0.0893 -0.230  Yes     
#> 8 Plaice Adult       0.0700  -0.178  -0.247  Yes

phi_pred_delta2 <- phi_pred_delta2 %>% 
  left_join(order %>% dplyr::select(Expected), by = "id2")
#> Adding missing grouping variables: `id2`
#> left_join: added one column (Expected)
#>            > rows only in x        0
#>            > rows only in y  (     0)
#>            > matched rows     16,000
#>            >                 ========
#>            > rows total       16,000

phi_pred_delta2$id2 <- factor(phi_pred_delta2$id2, levels = c(order$id2))
                                 
phi_pred_delta2 %>%
  ggplot(aes(bio_trend_delta, id2, fill = mean_clim, alpha = Expected)) +
  geom_vline(xintercept = 0, linetype = 2, color = "grey30") + 
  geom_density_ridges(color = NA) +
  scale_fill_manual(values = c("steelblue2", "tomato3")) +
  scale_alpha_manual(values = c(0.4, 0.8)) +
  labs(y = "", x = "\u0394 biotic velocity with a 1 st.dev\ndecline in \u03C6 velocity", fill = "Mean climate") +
  theme_plot(base_size = 16) +
  theme(legend.direction = "vertical") +
  NULL
#> Picking joint bandwidth of 0.0071


ggsave("figures/mi_velo_delta.pdf", width = 17, height = 17, units = "cm", device = cairo_pdf)
#> Picking joint bandwidth of 0.0071

# phi_pred %>% 
#   mutate(id2 = paste(life_stage, species)) %>% 
#   group_by(group) %>% 
#   mutate(mean_clim = ifelse(mean_phi_ct == max(mean_phi_ct), "High \u03C6", "Low \u03C6")) %>% 
#   ungroup() %>% 
#   ggplot(aes(phi_vel_sc, est, color = mean_clim, fill = mean_clim)) +
#   geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.3, color = NA) +
#   geom_line(alpha = 0.8) +
#   labs(y = "Predicted Biotic velocity (km/year)",
#        x = "\u03C6 velocity (km/year)") +
#   scale_fill_manual(values = c("steelblue2", "tomato3"), name = "") +
#   scale_color_manual(values = c("steelblue2", "tomato3"), name = "") +
#   facet_grid2(life_stage ~ species, scales = "free_y", independent = "y") +
#   theme(aspect.ratio = 1) + 
#   NULL

#ggsave("figures/mi_velo_conditional.pdf", width = 17, height = 17, units = "cm", device = cairo_pdf)

Fit sdmTMB models of biomass velocities as a function of temperature velocities, given mean temperature (ooxygen would be cool to but I think perhaps its too correlated with temperature, cor = 0.8 ish)

pred_velo_temp <- list()
pred_velo_temp2 <- list()

for(i in unique(d$group)){
  
  dd <- filter(d, group == i)
  
  mesh <- make_mesh(dd, c("X", "Y"), n_knots = 200)
  plot(mesh)
  
  m <- sdmTMB(
    bio_vel ~ temp_vel_sc * mean_temp_ct + mean_log_biomass_ct + quarter_f,
    mesh = mesh,
    data = dd,
    family = gaussian(), 
    spatial = "on",
    spatiotemporal = "off",
    control = sdmTMBcontrol(newton_loops = 1)
  )
  
  sanity(m) 
  
  # Predict
  low_mean_clim <- as.numeric(quantile(dd$mean_temp_ct, prob = 0.05))
  high_mean_clim <- as.numeric(quantile(dd$mean_temp_ct, prob = 0.95))
  mean_log_biomass_ct <- mean(dd$mean_log_biomass_ct)
  
  nd <- data.frame(expand.grid(temp_vel_sc = c(0, 1),
                               mean_temp_ct = seq(low_mean_clim,
                                                  high_mean_clim,
                                                  length.out = 2))) %>% 
    mutate(mean_log_biomass_ct = mean_log_biomass_ct,
           X = mean(dd$X),
           Y = mean(dd$Y),
           quarter_f = as.factor(4))
   
  nsim <- 1000
  pred <- predict(m, newdata = nd, nsim = nsim)

  pred <- data.frame(
    mean_temp_ct = rep(c(low_mean_clim, high_mean_clim), each = nsim),
    est_0 = c(pred[1, ], pred[3, ]),
    est_change = c(pred[2, ], pred[4, ])) # note different order here than phi

  pred <- pred %>% 
    mutate(delta_bio_vel = est_change - est_0, # if bio trend is positive, bio increases with declines in phi
           mean_clim = ifelse(mean_temp_ct == low_mean_clim, 
                              "low_temperature", "high_temperature"),
           group = i)
  
  pred_velo_temp[[i]] <- pred
  
  # Predict 2 (max and min warming, what's the expected?)
  # nd2 <- data.frame(expand.grid(temp_vel_sc = c(quantile(dd$temp_vel_sc, probs = 0.01),
  #                                               quantile(dd$temp_vel_sc, probs = 0.99)),
  #                               mean_temp_ct = seq(quantile(dd$mean_temp_ct, probs = 0.01),
  #                                                  quantile(dd$mean_temp_ct, probs = 0.99),
  #                                                  length.out = 2))) %>% 
  #   mutate(mean_log_biomass_ct = 0,
  #          X = mean(dd$X),
  #          Y = mean(dd$Y),
  #          quarter_f = as.factor(4))
  # 
  # pred2 <- predict(m, newdata = nd2) %>% mutate(group = i)
  # 
  # pred_velo_temp2[[i]] <- pred2

}
#> filter: removed 70,544 rows (80%), 17,341 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#>         new variable 'X' (double) with one unique value and 0% NA
#>         new variable 'Y' (double) with one unique value and 0% NA
#>         new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'delta_bio_vel' (double) with 2,000 unique values and 0% NA
#>         new variable 'mean_clim' (character) with 2 unique values and 0% NA
#>         new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 72,173 rows (82%), 15,712 rows remaining

#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#>         new variable 'X' (double) with one unique value and 0% NA
#>         new variable 'Y' (double) with one unique value and 0% NA
#>         new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'delta_bio_vel' (double) with 2,000 unique values and 0% NA
#>         new variable 'mean_clim' (character) with 2 unique values and 0% NA
#>         new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 86,033 rows (98%), 1,852 rows remaining

#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#>         new variable 'X' (double) with one unique value and 0% NA
#>         new variable 'Y' (double) with one unique value and 0% NA
#>         new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'delta_bio_vel' (double) with 2,000 unique values and 0% NA
#>         new variable 'mean_clim' (character) with 2 unique values and 0% NA
#>         new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 86,343 rows (98%), 1,542 rows remaining

#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#>         new variable 'X' (double) with one unique value and 0% NA
#>         new variable 'Y' (double) with one unique value and 0% NA
#>         new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'delta_bio_vel' (double) with 2,000 unique values and 0% NA
#>         new variable 'mean_clim' (character) with 2 unique values and 0% NA
#>         new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 64,646 rows (74%), 23,239 rows remaining

#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#>         new variable 'X' (double) with one unique value and 0% NA
#>         new variable 'Y' (double) with one unique value and 0% NA
#>         new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'delta_bio_vel' (double) with 2,000 unique values and 0% NA
#>         new variable 'mean_clim' (character) with 2 unique values and 0% NA
#>         new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 79,901 rows (91%), 7,984 rows remaining

#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✖ `b_j` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> ✖ `b_j` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> ✖ `b_j` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#>         new variable 'X' (double) with one unique value and 0% NA
#>         new variable 'Y' (double) with one unique value and 0% NA
#>         new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'delta_bio_vel' (double) with 2,000 unique values and 0% NA
#>         new variable 'mean_clim' (character) with 2 unique values and 0% NA
#>         new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 79,755 rows (91%), 8,130 rows remaining

#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#>         new variable 'X' (double) with one unique value and 0% NA
#>         new variable 'Y' (double) with one unique value and 0% NA
#>         new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'delta_bio_vel' (double) with 2,000 unique values and 0% NA
#>         new variable 'mean_clim' (character) with 2 unique values and 0% NA
#>         new variable 'group' (character) with one unique value and 0% NA
#> filter: removed 75,800 rows (86%), 12,085 rows remaining

#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'mean_log_biomass_ct' (double) with one unique value and 0% NA
#>         new variable 'X' (double) with one unique value and 0% NA
#>         new variable 'Y' (double) with one unique value and 0% NA
#>         new variable 'quarter_f' (factor) with one unique value and 0% NA
#> mutate: new variable 'delta_bio_vel' (double) with 2,000 unique values and 0% NA
#>         new variable 'mean_clim' (character) with 2 unique values and 0% NA
#>         new variable 'group' (character) with one unique value and 0% NA


pred_velo_temp <- dplyr::bind_rows(pred_velo_temp) %>% 
  separate(group, c("species", "life_stage"),
           sep = "_", remove = FALSE) %>% 
  mutate(species = str_to_title(species),
         life_stage = str_to_title(life_stage))
#> mutate: changed 16,000 values (100%) of 'species' (0 new NA)
#>         changed 16,000 values (100%) of 'life_stage' (0 new NA)

# pred_velo_temp2 <- dplyr::bind_rows(pred_velo_temp2) %>% 
#   separate(group, c("species", "life_stage"),
#            sep = "_", remove = FALSE) %>% 
#   mutate(species = str_to_title(species),
#          life_stage = str_to_title(life_stage))

write_csv(pred_velo_temp, "output/pred_temp_velo.csv")
pred_velo_temp2 <- pred_velo_temp %>%
  mutate(mean_clim = ifelse(mean_clim == "low_temperature", "Cold", "Warm")) %>%
  mutate(id2 = paste(species, life_stage))
#> mutate: changed 16,000 values (100%) of 'mean_clim' (0 new NA)
#> mutate: new variable 'id2' (character) with 8 unique values and 0% NA
  
order <- pred_velo_temp2 %>%
  group_by(id2, mean_clim) %>% 
  summarise(mean_delta = mean(delta_bio_vel)) %>% 
  pivot_wider(names_from = mean_clim, values_from = mean_delta) %>% 
  mutate(diff = `Warm` - `Cold`) %>% 
  arrange(desc(diff)) %>% 
  mutate(Expected = ifelse(Cold > Warm, "Yes", "No"))
#> group_by: 2 grouping variables (id2, mean_clim)
#> summarise: now 16 rows and 3 columns, one group variable remaining (id2)
#> pivot_wider: reorganized (mean_clim, mean_delta) into (Cold, Warm) [was 16x3, now 8x3]
#> mutate (grouped): new variable 'diff' (double) with 8 unique values and 0% NA
#> mutate (grouped): new variable 'Expected' (character) with 2 unique values and 0% NA

order
#> # A tibble: 8 × 5
#> # Groups:   id2 [8]
#>   id2                  Cold     Warm     diff Expected
#>   <chr>               <dbl>    <dbl>    <dbl> <chr>   
#> 1 Flounder Adult    -0.0226  0.106    0.128   No      
#> 2 Dab Adult          0.0289  0.0553   0.0264  No      
#> 3 Dab Juvenile       0.0142  0.0313   0.0171  No      
#> 4 Plaice Juvenile    0.0424  0.0509   0.00852 No      
#> 5 Plaice Adult      -0.0250 -0.0382  -0.0132  Yes     
#> 6 Cod Adult          0.0330 -0.00467 -0.0377  Yes     
#> 7 Flounder Juvenile  0.0488 -0.0242  -0.0730  Yes     
#> 8 Cod Juvenile       0.0410 -0.0358  -0.0769  Yes

pred_velo_temp2 <- pred_velo_temp2 %>% 
  left_join(order %>% dplyr::select(Expected), by = "id2")
#> Adding missing grouping variables: `id2`
#> left_join: added one column (Expected)
#>            > rows only in x        0
#>            > rows only in y  (     0)
#>            > matched rows     16,000
#>            >                 ========
#>            > rows total       16,000

pred_velo_temp2$id2 <- factor(pred_velo_temp2$id2, levels = c(order$id2))
                                 
pred_velo_temp2 %>%
  ggplot(aes(delta_bio_vel, id2, fill = mean_clim, alpha = Expected)) +
  geom_vline(xintercept = 0, linetype = 2, color = "grey30") + 
  geom_density_ridges(color = NA) +
  scale_fill_manual(values = c("steelblue2", "tomato3")) +
  scale_alpha_manual(values = c(0.4, 0.8)) +
  labs(y = "", x = "\u0394 biotic velocity with a 1 st.dev\nincrease in temperature velocity", fill = "Mean temperature") +
  theme_plot(base_size = 16) +
  theme(legend.direction = "vertical") +
  NULL
#> Picking joint bandwidth of 0.00598


ggsave("figures/temp_velo_delta.pdf", width = 15, height = 17, units = "cm", device = cairo_pdf)
#> Picking joint bandwidth of 0.00598
knitr::knit_exit
#> function (append, fully = TRUE) 
#> {
#>     if (missing(append)) 
#>         append = if (out_format(c("latex", "sweave", "listings"))) 
#>             "\\end{document}"
#>         else if (out_format("html")) 
#>             "</body>\n</html>"
#>         else ""
#>     .knitEnv$terminate = append
#>     .knitEnv$terminate_fully = fully
#>     invisible()
#> }
#> <bytecode: 0x7fa2ad793708>
#> <environment: namespace:knitr>
# pred_velo_temp2 %>%
#   mutate(group2 = paste(species, life_stage, sep = "")) %>% 
#   group_by(group) %>% 
#   mutate(mean_clim = ifelse(mean_temp_ct == min(mean_temp_ct), "Low temperature", "High temperature"),
#          warming = ifelse(temp_vel_sc == max(temp_vel_sc), "Max warming", "Min warming")) %>%
#   ungroup() %>% 
#   ggplot(aes(est, group, color = mean_clim)) +
#   geom_point(alpha = 0.7, size = 3) +
#   scale_color_manual(values = c("red", "blue")) +
#   geom_vline(xintercept = 0, linetype = 2) +
#   facet_wrap(~factor(warming, levels = c("Min warming", "Max warming"))) +
#   scale_fill_gradient2()
# 
# # Filter only max warming
# pred_velo_temp2 %>%
#   mutate(group2 = paste(species, life_stage, sep = "")) %>% 
#   group_by(group) %>% 
#   mutate(mean_clim = ifelse(mean_temp_ct == min(mean_temp_ct), "Low temperature", "High temperature"),
#          warming = ifelse(temp_vel_sc == max(temp_vel_sc), "Max warming", "Min warming")) %>%
#   ungroup() %>% 
#   filter(warming == "Max warming") %>% 
#   ggplot(aes(est, group, color = mean_clim)) +
#   geom_point(alpha = 0.7, size = 3) +
#   scale_color_manual(values = c("red", "blue")) +
#   geom_vline(xintercept = 0, linetype = 2) +
#   facet_wrap(~factor(warming, levels = c("Min warming", "Max warming"))) +
#   scale_fill_gradient2()